Format Cetaceans Data

Creation - jeremy.andreoletti@ens.fr - 13/04/2020

Load occurrence dataset

cannot open file 'data-cetaceans/cetacea_pbdb_11_May_2020.csv': No such file or directoryError in file(file, "rt") : cannot open the connection

Reorder accepted ranks according to classification standard.

 [1] "family"         "genus"          "infraorder"     "order"          "species"        "subfamily"      "suborder"      
 [8] "subspecies"     "superfamily"    "unranked clade"

Load the list of occurrences with morphological information

                      Taxon    composite..n.y.     occ.data.based.on        Specimen  pbdb_specimen_no pbdb_occurence.number
 Aetiocetus_cotylalveus  : 1   Min.   :0.0000               :60      AMNH 9485  : 1   Min.   : 25492   Min.   :  68135      
 Agorophius_pygmaeus     : 1   1st Qu.:0.0000   CMM-V-15    : 1      CASG 66660 : 1   1st Qu.: 25706   1st Qu.: 461107      
 Albertocetus_meffordorum: 1   Median :0.0000   MGGC 8548   : 1      CCNHM-101  : 1   Median : 25922   Median : 487310      
 Albireo_whistleri       : 1   Mean   :0.3768   MNHN SAS 933: 1      ChM PV4256 : 1   Mean   : 35167   Mean   : 631596      
 Archaeodelphis_patrius  : 1   3rd Qu.:1.0000   UCMP 83790  : 1      ChM PV4844 : 1   3rd Qu.: 26035   3rd Qu.: 763082      
 Ashleycetus_planicapitis: 1   Max.   :1.0000   USNM 10484  : 1      CHM_PV_4253: 1   Max.   :146940   Max.   :1360382      
 (Other)                 :63                    (Other)     : 4      (Other)    :63   NA's   :56                            
                                                                                                                                                                             notes   
                                                                                                                                                                                :62  
 actually this info relates to the holotype, which is MGGC 8608, but it seems that MGGC 8599 has not been entered                                                               : 1  
 It is unclear which specimen is in the pbdb, both are from the Pisco Formation. This one is the referred specimen, not the holotype (but ir is more complete than the holotype): 1  
 two specimen numbers in pbdb, but both relate to same specimen                                                                                                                 : 1  
 two specimens in database, but one comes from USA, so I selected the one from Peru based on collection info                                                                    : 1  
 two specimens in pbdb, unlcear which is which because they belong to the same horizon                                                                                          : 1  
 (Other)                                                                                                                                                                        : 2  

Combine those taxa with extant taxa to get all the species included in the tree.

                        taxon         min         max   
 Balaena_mysticetus        : 1   Min.   :0   Min.   :0  
 Balaenoptera_acutorostrata: 1   1st Qu.:0   1st Qu.:0  
 Balaenoptera_bonaerensis  : 1   Median :0   Median :0  
 Balaenoptera_borealis     : 1   Mean   :0   Mean   :0  
 Balaenoptera_brydei       : 1   3rd Qu.:0   3rd Qu.:0  
 Balaenoptera_edeni        : 1   Max.   :0   Max.   :0  
 (Other)                   :82                          

                      taxon         min              max       
 Aetiocetus_cotylalveus  : 1   Min.   : 2.588   Min.   : 3.60  
 Agorophius_pygmaeus     : 1   1st Qu.: 7.246   1st Qu.:11.62  
 Albertocetus_meffordorum: 1   Median :13.820   Median :15.97  
 Albireo_whistleri       : 1   Mean   :14.357   Mean   :18.28  
 Archaeodelphis_patrius  : 1   3rd Qu.:23.030   3rd Qu.:28.10  
 Ashleycetus_planicapitis: 1   Max.   :41.300   Max.   :47.80  
 (Other)                 :63                                   

                        taxon          min             max        
 Balaena_mysticetus        :  1   Min.   : 0.00   Min.   : 0.000  
 Balaenoptera_acutorostrata:  1   1st Qu.: 0.00   1st Qu.: 0.000  
 Balaenoptera_bonaerensis  :  1   Median : 0.00   Median : 0.000  
 Balaenoptera_borealis     :  1   Mean   : 6.31   Mean   : 8.035  
 Balaenoptera_brydei       :  1   3rd Qu.:11.62   3rd Qu.:13.820  
 Balaenoptera_edeni        :  1   Max.   :41.30   Max.   :47.800  
 (Other)                   :151                                   

Check that the names in the datasets are exactly the same as in the list above.

[1] TRUE

Idem at the generic level for the genus analysis.

                        taxon         min         max   
 Balaena_mysticetus        : 1   Min.   :0   Min.   :0  
 Balaenoptera_physalus     : 1   1st Qu.:0   1st Qu.:0  
 Berardius_bairdii         : 1   Median :0   Median :0  
 Caperea_marginata         : 1   Mean   :0   Mean   :0  
 Cephalorhynchus_heavisidii: 1   3rd Qu.:0   3rd Qu.:0  
 Delphinapterus_leucas     : 1   Max.   :0   Max.   :0  
 (Other)                   :35                          
                      taxon         min              max       
 Aetiocetus_cotylalveus  : 1   Min.   : 2.588   Min.   : 3.60  
 Agorophius_pygmaeus     : 1   1st Qu.: 7.246   1st Qu.:11.62  
 Albertocetus_meffordorum: 1   Median :13.820   Median :15.97  
 Albireo_whistleri       : 1   Mean   :14.960   Mean   :18.95  
 Archaeodelphis_patrius  : 1   3rd Qu.:23.030   3rd Qu.:28.10  
 Ashleycetus_planicapitis: 1   Max.   :41.300   Max.   :47.80  
 (Other)                 :56                                   
                        taxon         min              max        
 Balaena_mysticetus        : 1   Min.   : 0.000   Min.   : 0.000  
 Balaenoptera_physalus     : 1   1st Qu.: 0.000   1st Qu.: 0.000  
 Berardius_bairdii         : 1   Median : 5.333   Median : 7.246  
 Caperea_marginata         : 1   Mean   : 9.005   Mean   :11.409  
 Cephalorhynchus_heavisidii: 1   3rd Qu.:15.970   3rd Qu.:20.440  
 Delphinapterus_leucas     : 1   Max.   :41.300   Max.   :47.800  
 (Other)                   :97                                    

Check that the names in the datasets are exactly the same as in the list above.

[1] TRUE

Remove those occurrences from our initial dataset to avoid redundancy.

Explore the dataset

Repartition through time

Full fossil record

\(\to\) Numerous occurrences seem to have the same age interval so in order to avoid clusters let’s draw them uniformly in their stratigraphic range rather than taking the mean.

\(\to\) The repartition seems much smoother now.

Subsampling

These occurrences are too numerous for our current implementation, let’s subsample a fraction of them for now.

\(\to\) The distribution looks similar, with some noise due to higher variance with smaller sample.

Repartition among accepted ranks

Pie chart

\(\to\) Half of the occurrences are identified at the level of the secies and 1/3 at the genus or family.

Some clades are unranked :

  Chaeomysticeti          Neoceti Panphyseteroidea       Pelagiceti  Platanidelphidi       Squaloceti 
              26                2                1                1                1                1 

Time repartition by rank

\(\to\) Apparent huge cluster of occurrences in recent times, with very precise dating = Artefact due to the “Pull of the Recent” effect ? \(\to\) We decided to remove all Late Pleistocene and Holocene occurrences (thus setting the \(\omega\)-sampling to 0) in order to avoid this bias.

\(\to\) We observe similar trends at each rank, with peaks at ~15My and ~5My.

Redundancy of occurrences with the same accepted name

\(\to\) ~Half of species/genera/subfamilies have only one specimen by accepted name, but it could go up to ~50 within the same species and ~200 occurrences within the same suborder. Those differences will have to be corrected because in our model all species are upposed to have the same abundance (identical sampling rates among branches).

\(\implies\) Our goal now will be to correct this abundance bias.

Time intervals = stratigraphic age uncertainty

Minimum and maximum stratigraphic limits

\(\to\) Most species have a early but not a late stratigraphic limit.

Minimum and maximum ages

Time ranges = duration of the time intervals

Count occurrences by accepted name

\(\to\) Some occurrences have too much age uncertainty, they risk to artificially increase species durations.

Remove occurrences with highly uncertain dating (range > 10My)

[1] 3505  117

Most of occurrences show less than 10 My age uncertainty, let’s try to keep only these ones.

\(\to\) The removal of highly uncertain occurrences seems to be only a little biased, even if uncertainty globally increases with age.

\(\to\) Some species (or other ranks) have several occurrences with several time ranges, let’s combine them into a unique range covering all the others.

Combined time ranges = unique time range for occurrences with the same name (without the biggest ones)

Occurrence density

Density distributions

\(\to\) Density logically increases as taxa ranks increase.

Correlation between time range and age

If we want to correct species abundance differences based on the number of occurrences in the time range (“density”), those factors should not depend on time in order to avoid penalizing periods with higher densities.

\(\to\) It seems that age range varies importantly with time, but when taking the full combined range into account the correlation seems quite weak after the first million years.

Let’s look at the density directly, because this is what is interesting us directly.

\(\to\) The density based on combined ranges is much less time-dependant than the density based on initial range ages. We will therefore use the combined density for our corrections.

Sub-sampling of occurrences with a normalized density along the combined ranges

Compare densities for single vs. combined ranges.

\(\to\) Densities are smaller and more concentrated with the combined ranges (larger time span + less ranges in total because of the collepses into unique ranges).

Compare densities by accepted name count (species only)

Let’s focus now on the occurrences accepted at the species level because they are the one for which we can correct the abundance bias by subsampling the most concentrated combined intervals.

\(\to\) Their is a huge span of densities driven by the number of occurrences for the same species that we can reduce by subsampling the most concentrated intervals.

Impact of correcting subsampling on density distributions (species only)

\(\to\) Subsampling successfully reduces the density span from 2 to 1 order of magnitude.

Impact of subsampling on occurrences repartition (species only)

See what our distributions loook like after subsampling :

\(\to\) Some highly dense cluster became much more similar to the others.

If we superpose these 2 plots :

\(\to\) We get the new species occurrence repartition after subsampling correction, that could be used for doing inference with the occurrence birth-death model.

New developments

Compare with a Poisson sampling process

In order to check if the data fit our assumptions of constant-fossilisation-rate Poisson sampling we compare the observed occurrences distributions with the expected ones. Specifically, we will look at the number of taxa represented by 1, 2, 3, … occurrences and the one that we would expect for a given distribution of combined age ranges (as a proxy for species duration).

In a Poisson process with occurrence sampling rate \(\omega\) and for a given time interval of length \(t\), the probability of observing \(N_t = k\) occurrences is given by the Poisson distribution of mean of parameter \(\omega\times t\) :

\[\mathbb{P}(N_t=k) = e^{-\omega t}\frac{(\omega t)^k}{k!}\]

So in order to have the absolute probability of observing \(N_0 = n\) occurrences we have to integrate over the full distribution of age ranges \(t\), called \(f(t)\) :

\[\mathbb{P}(N_0 = n) = \int_t P(N_t = n) f(t) dt = \int_t e^{-\omega t}\frac{(\omega t)^n}{n!} f(t) dt\]

First, approximate this distribution :

Then integrate and plot the expected distribution for a given omega :

Finally, try to find an \(\omega\) value that approximately fits the first points (the least affected by oversampling biases) and check if the other points follow the expected curve :

\(\implies\) Initial observations really do not fit the expectations, species with more than 5 occurrences must remain very rare ! But our subsampling seems to correct most of this bias.

However, this method requires to make several arbitrary choices that may introduce new biases so we will instead subsample at other levels (palaeontological collection, geological formation). In each case only one occurrence will be sampled for the similarly identified, a process we will refer to as aggregating these occurrences according to the chosen factor.

Aggregate similarly identified occurrences in each collection

\(\to\) Each collection corresponds to a unique time interval (inferred from the unique age mean).

In order to reduce the abundance bias, we may keep only one occurrence for each collection :

\(\to\) Not enough occurrences are removed to make a sufficient difference. If we look at the collection with the highest number of occurrences :

        Aprixokogia kelloggi                      Balaena   Balaenoptera acutorostrata              Balaenopteridae 
                           1                            2                            2                            2 
                   Balaenula      Bohaskaia monodontoides                Cetotheriinae               Delphinapterus 
                           2                            1                            1                            3 
                 Delphinidae                    Delphinus                 Globicephala          Gricetoides aurorae 
                           2                            2                            2                            1 
              Herpetocetinae       Herpetocetus sendaicus Herpetocetus transatlanticus              Kogia breviceps 
                           1                            1                            1                            1 
                    Kogiidae                     Kogiinae          Kogiopsis floridana               Lagenorhynchus 
                           1                            2                            1                            2 
    Lagenorhynchus harmatuki                    Megaptera      Mesoplodon longirostris                      Monodon 
                           1                            2                            2                            1 
    Ninoziphius platyrostris                Orycterocetus       Physeter macrocephalus                 Physeteridae 
                           3                            1                            1                            1 
                Physeterinae           Physeterula dubusi                  Plesiocetus        Pliopontos littoralis 
                           2                            1                            1                            1 
                  Pontoporia                Pontoporiidae                    Pseudorca                  Scaldicetus 
                           1                            2                            2                            1 
                    Stenella                Stenella rayi                     Tursiops          Ziphius cavirostris 
                           2                            1                            2                            2 

\(\to\) There are very few redundancies among the accepted names in collections so aggregating those won’t reduce the abundance bias.

Instead, we may try to aggregate occurrences with the same accepted name at the level of the geological formation (ie subsample only one for each).

Aggregate similarly identified occurrences in each formation

\(\to\) For most formations each species seems to be restricted to only one age, so as expected we are not loosing too much information when aggregating them to a single occurrence.

\(\to\) In that case the sub-sampling is big enough to hope correcting our bias.

\(\to\) The bias can be corridered as completely corrected.

Removed 2 rows containing missing values (geom_bar).

Removed 2 rows containing missing values (geom_bar).

Removed 2 rows containing missing values (geom_bar).Removed 2 rows containing missing values (geom_bar).

\(\to\) Comparing with the initial occurrences distribution and with the distribution after our first sub-sampling it appears that both methods lead to very similar distributions. This conforts us about the robustness of those approaches.

However, we have to take into account the occurrences that do not have any indicated geological formation to subsample them separately. To approximate geological formation we chose to proceed to the aggregation based on the combination of the country and the early stratigraphic interval.

Aggregate occurrences without formation by country + early interval

More occurrences remain after aggregating with this new method. Let’s compare again with the theoretical distribution :

\(\to\) The correspondance is still good, except for two taxa :

   Scaldicetus grandis Schizodelphis sulcatus 
                    14                     16 

Removed 2 rows containing missing values (geom_bar).

Removed 2 rows containing missing values (geom_bar).

Removed 2 rows containing missing values (geom_bar).Removed 2 rows containing missing values (geom_bar).

\(\to\) Comparing with the initial occurrences distribution and with the distribution after our first sub-sampling it appears that both methods lead to very similar distributions. This conforts us about the robustness of those approaches.

Replacing coutries by gelogical plates seems to make more sens from a palaeontological perspective, so let’s try it.

Aggregate occurrences without formation by geoplate + early interval

More occurrences remain after aggregating with this new method. Let’s compare again with the theoretical distribution :

\(\to\) The correspondance is still good, except for two taxa :

   Scaldicetus grandis Schizodelphis sulcatus 
                    13                     15 

Removed 2 rows containing missing values (geom_bar).

Removed 2 rows containing missing values (geom_bar).

Removed 2 rows containing missing values (geom_bar).Removed 2 rows containing missing values (geom_bar).

\(\to\) Delimiting by geological plates (+ age) instead of countries (+ age) leads to similar distributions, so we will keep it.

Check that the sampling methods do not introduce biases in the repartition between Odontoceti and Mysticeti

Mystecetes are usually larger than odontocetes, and size is associated with a wider geographic range so since we are subsampling occurrences according to geological formation we may be biasing our data towards more widespread species, therefore towards mystecetes.

Look at the families first :

\(\to\) The family repartitions vary a bit after subsampling, but because of the limited number of species by family the fact that we corrected the oversampling of some species could have a disproportionate effect.

Let’s look rather at a higher taxonomic rank, by importing the topology of cetacean families (from Marx et al. 2016) :

`data_frame()` is deprecated as of tibble 1.1.0.
Please use `tibble()` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.

\(\to\) There is a smaller proportion of archeocete occurrences after either sampling, because a huge cluster is subsampled around 35 My ago, but this effect is expected. However, for the Mysticeti vs. Odotoceti there is no huge apparent bias, especially with the aggregating method.

Faster genus-level analysis

More occurrences remain after aggregating with this new method. Let’s compare again with the theoretical distribution :

\(\to\) The correspondance is much less improved than with species aggregation because some genera have a lot of occurrences due to their high number of species species :

Balaenoptera   Kentriodon   Mesoplodon  Scaldicetus    Squalodon 
          53           27           30           41           40 

Removed 2 rows containing missing values (geom_bar).

Removed 2 rows containing missing values (geom_bar).

Removed 2 rows containing missing values (geom_bar).Removed 2 rows containing missing values (geom_bar).

\(\to\) Delimiting by geological plates (+ age) instead of countries (+ age) leads to similar distributions, with a bit more occurrences, so we will keep it.

Conclusions

Achievements :

• It seems possible to adequately reduce the abundance bias by subsampling the most concentrated intervals \(\to\) species only
• Using combined ranges by species appears to be more robust \(\to\) to be confirmed
• Very recent samples may have been dated with a more precise method and contain much more fossils, so they should be removed or treated separatadely \(\to\) additional information needed

Open questions :

• What about other accepted ranks ?
  1. The problem is that differences in the number of occurrences at higher ranks could be due to differences in individual abundances inside species or due to differences in the number of species inside that group.
  2. A solution could be to look at the number of species by group based on the indicated species, and include it in the bias correction : homogeneizing the number of occurrences / time unit / number sf species in the group \(\to\) additional data required (ranks classification)
• Why do most occurrences miss a late stratigraphic limit ?
• Some occurrences have very huge time intervals \(\to\) Was is a good idea to remove those >10My ? Should we remove more of them (>5My) ?